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ABSTRACT 

Context. During the formation of a star, material is ejected along powerful jets that impact the ambient material. This outflow regulates 
star formation by e.g. inducing turbulence and heating the surrounding gas. Understanding the associated shocks is therefore essential 
to the study of star formation. 

Aims. We present comparisons of shock models with CO, H 2 , and SiO observations in a ‘pure’ shock position in the BHR71 bipolar 
outflow. These comparisons provide an insight into the shock and pre-shock characteristics, and allow us to understand the energetic 
and chemical feedback of star formation on Galactic scales. 

Methods. New CO (Ap = 16, 11, 7, 6,4, 3) observations from the shocked regions with the SOFIA and APEX telescopes are presented 
and combined with earlier H 2 and SiO data (from the Spitzer and APEX telescopes). The integrated intensities are compared to a grid 
of models that were obtained from a magneto-hydrodynamical shock code which calculates the dynamical and chemical stmcture of 
these regions combined with a radiative transfer module based on the ‘large velocity gradient’ approximation. 

Results. The CO emission leads us to update the conclusions of our previous shock analysis: pre-shock densities of 10'' cm“^ and 
shock velocities around 20-25 km s“* are still constrained, but older ages are inferred (~4000 years). 

Conclusions. We evaluate the contribution of shocks to the excitation of CO around forming stars. The SiO observations are com¬ 
patible with a scenario where less than 4% of the pre-shock SiO belongs to the grain mantles. We infer outflow parameters: a mass 
of 1.8 X 10“^ Mq was measured in our beam, in which a momentum of 0.4 Mq km s“' is dissipated, corresponding to an energy of 
4.2 X 10''^ erg. We analyse the energetics of the outflow species by species. Comparing our results with previous studies highlights 
their dependence on the method: H 2 observations only are not sufficient to evaluate the mass of outflows. 

Key words. Stars: formation - ISM: jets and outflows - ISM: individual objects: BHR71 - Submillimeter: ISM - Infrared: ISM 


1. Introduction 

Molecular shocks are ubiquitous in the interstellar medium 
(ISM) of our Galaxy. They can be associated with the forma¬ 
tion of ‘ridges’ at the convergence region of molecular clouds 
(e.g. W43, Nguyen-Lu’o’ng et al. 2013| l, to the jets and out¬ 
flow systems related to the formation of low- to high-mass stars 
(from low, e.g. |Gomez-Ruiz et al.||2013| to high, |Leurini et ah] 
|2013 through intermediate, Gomez-Ruiz et al.|2()12 for exam¬ 
ples; also see |Arce et al.||2007| [Frank et al.||2014| or jTan et al.| 
|2014| for reviews), or to supernova remnants (SNRs) interacting 
with molecular clouds (e.g. W44, |Anderl et al.|2014 l. However, 
an analysis of these shocks is challenging because several phys¬ 
ical processes are contributing to the observed emission. The 
W43 ridges are illuminated by an energetic UV radiation field 
coming from neighbouring Hit regions of star clusters ( Bally] 
et al.||2010l. Studying young stellar objects (YSO) by pointing 


jet al.|2013] l, infall processes, and UV illumination (the latter be¬ 
ing all the more important when the mass of the forming object 
I large, e.g. Visser et al. 2012 San Jose-Garcla et al.|2013|l. Out¬ 


flows from massive star-forming regions are also illuminated by 


strong radiation in the X-ray regime (e.g. W28 A2, Rowell et al. 
2010|l. Similarly, old SNRs are also often subject to X-rays, and 


also to y-ray emission which are the signature of the acceleration 
of particles that locally took place shortly after the supernova ex¬ 
plosion (e.g. Hewitt et al.|2009]l. 


‘Pure’ molecular shock regions can be defined as regions 
where the physics and chemistry are dominantly driven by 
shocks. Studying these regions, therefore, is of crucial impor¬ 
tance to investigating the feedback they exert on their envi¬ 
ronment, whether this feedback is energetic or chemical. From 
the energetic point of view, these studies allow us to assess 
the contribution of shocks to the excitation of e.g. CO on 


on the central protostar leads in practice to a study of the com¬ 
bination of ejection shocks (that can be multiple, e.g. jKristens^ 


galactic scales, as observed by Herschel (in NGC1068: Hailey- j 
jPunsheath et al.|2012| M82: |Kamenetzky et al.|20T2l NGC6240 
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and Mrk231: |Meijerink et alT 2013 or NGC253: [Rosenberg et aL] 
2014a|and Arp299: Rosenberg et al.||2014bl from the inside of 


the Galaxy. From the chemical point of view, these studies al¬ 
low us to investigate the formation paths of water (e.g. [Leurinij 
et al.|2014bll or more complex molecules (e.g. complex organic 


Belloche et al.|20r3) l, and to understand their presence in 


ones, 

planetary systems, for instance. These ‘pure’ shock regions are 
not numerous. The most remarkable and well-studied pure shock 
region is the B1 knot of the LI 157 bipolar outflow. Sufficiently 
distant from the central protostar, this region remains uncontam¬ 
inated by infall or irradiation processes, and has been the sub¬ 
ject of a number of studies dedicated to studying its energetics 
or chemical composition (Gusdorf et al.||2008b[ Codella et ak] 


2010[|Fl ower & Pineau d es Forets|2012[|Benedettini et al.|2013t 

Busquet et al.|2014[]Podio et al.|2014| l. In particular, LI 157 was 


mapped by the GREAT (german receiver for astronomy at ter¬ 
ahertz frequencies) receiver onboard SOFIA (Stratospheric Ob¬ 
servatory for Infrared Astronomy) in CO (11-10), but the low 
signal-to-noise ratio prevented jEisloffel et al.|p012| l from a thor¬ 
ough study of its energetics or chemistry. This article focus on 
the analysis of a similar shock position in the BHR71 outflow, 
the ‘southern twin’ of LI 157. Because of its southern location, 
BHR71 is an ideal target to be observed with ALMA (Ata¬ 
cama Large Millimeter/submillimeter Array), which will never 
be done for LI 157 because of its high northern location. This ar¬ 
ticle is organised as follows: sectionj^presents our observations. 
Sectionj^presents the new CO data we made use of in our analy¬ 
sis, while the existing H 2 and SiO data is described in section]^ 
The results of shock modelling and their further use is exposed 
in section]^ and section [^summarises our findings. 


2. The BHR71 bipolar outflow 

2.1. Previous work 


The double bipolar outflow BHR71 (jBourke et al.|1997| Myers 


|& Mardones|199^ jParise et al.|2006| l lies close to the plane of 

the sky. It is powered by two sources IRS 1 and IRS2, separated 
by ~3400 AU ( Bourke et al.||200T] l, of luminosities 13.5 and 


et al. 

2011 

et al. 

2011 


hereafter N09 and Gial 1), and H 2 plus SiO ( [Gusdorf 


0.5 Lq ([Chen et al.[|2008[ l, relatively nearby (~200 pc, [Bourke [ 
[et al.|1995a|l. The Resent work builds on previous studies of the 


hereafter G11). In their work, N09 mostly described 
the Spitzer observations of the outflow. The InfraRed Spectro¬ 
graph (IRS) was used to map the inner part of the outflow in the 
pure rotational transitions of H 2 as well as in Fe ii and S i transi¬ 
tions. A region corresponding to approximately half the length of 
the outflow was covered by these observations around the driving 
protostars. The results were compared to previous observations 
of the entire region by the InfraRed Array Camera (IRAC) on¬ 
board the same telescope, showing that 30% and 100% of the 
luminosity of bands 3 and 4 could be accounted for by H 2 lines 
emission, similar to LI 157. The pure rotational H 2 luminosity of 
the flow was estimated to be 4.4x10“^ Lq, less than 1/3 of that 
of LI 157, but measured only from the fraction of the BHR71 
outflow that was mapped. The H 2 mass above 100 K was con¬ 
strained to be around 2.5x10“^ Mq, 20 times less than in LI 157. 
However the H 2 densities for both outflows were constrained to 
~10^ * cm“^. In Giall, these IRS observations of H 2 were de¬ 
tailed line by line. An average column density of H 2 of around 
10^° cm“^ was extracted from the map. Two temperature com¬ 
ponents were apparent, at ~300 and ~1500 K. A non-local ther¬ 
modynamical equilibrium analysis was performed and yielded a 


high H 2 average density of a few 
nosity of the flow was estimated to be twice as high as the pure 
rotational lines, also based on previous observations of the rovi- 
brational lines emission presented in [Giannini et al.[ ( [2004| l. In a 
few positions where this was possible, the observations of pure 
rotational H 2 lines were combined with those of rovibrational 
lines to generate a more complete excitation diagram that was 
also compared to shock models, confirming the high value of H 2 
density. 

In Gil, the authors presented a shock-model analysis of 
various positions in the northern lobe of the BHR71 outflow. 
They used Spitzer observations of H 2 , and APEX observations 
of SiO to constrain shock models. Their analysis was hampered 
by the rather high beam-size of their SiO observations. How¬ 
ever, they identified two positions where a shock analysis was 
favourable: the ‘SiO peak’ and ‘SiO knot’. These positions lie at 
the apex of the inner bipolar structure of the outflow, relatively 
un-contaminated by envelope infall. They are far away from the 
protostars, and separated from them by the inner outflow cavity, 
hence as shielded as possible from their potential UV radiation. 
For all these reasons they are reminiscent of the L1157-B1 po¬ 
sition, and are considered in the present study as ‘pure shock’ 
positions. 

Figure [T] shows the entire BHR71 outflow as seen through 
these various datasets and highlights the particular positions as 
first introduced by |Giannini et al.| ( 2004| l: the so-called knots 5,6, 
and 8 are thus indicated, as well as the two Herbig-Haro objects 
HH320A/B and HH321A/B. The driving protostars IRS 1 and 
IRS 2 are also marked. The two extra positions used in G11, ‘SiO 
peak’ (of emission) and neighbouring ‘SiO knot’ are also shown 
at the apex of the upper inner lobe of the outflow. In this fig¬ 
ure, the previous and most recent observations can be compared 
(APEX, see section jL^ : intensity in the CO (6-5) (colours) and 
(3-2) (white contours) transitions integrated between -50 and 50 
km s"' in the left panel, and Spitzer IRAC (8 pm, colours) and 
IRS (H 2 0-0 S(5), white contours) maps overlaid with SiO (5—4) 
map (red contours in the upper lobe) in the right panel. The white 
contours then show how much of the outflow was observed by 
the IRS onboard Spitzer. 


2.2. APEX observations of CO 

APEX0 observations of the entire BHR71 outflow were con¬ 
ducted in 2013 and 2014. The analysis of the whole maps is 
out of the scope of the present work, and will be the subject 
of a forthcoming publication. In the present study, we used in¬ 
formation inferred from the full maps in the '^CO (3-2), '^CO 
(3-2), (4-3), (6-5), and (7-6) transitions, which we briefly de¬ 
scribe. The APEX observations towards BHR71 were conducted 
on several days: June 3 and 4, 2013, and June 28, 2014. We 
used the heterodyne receivers FLASH345, FLASH460 (First 
Light APE X Submillimeter Heterodyne receiver, [Heyminck[ 
et al. 2006 1 , and CHAMP^ (Carbon Heterodyne Array of the 
MPIfRT TKasemann et al.|[2006[ [Glisten et al.[[200'8l ), in com¬ 
bination with the MPIfR fast Eourier transform spectrometer 
backend (XEETS, [Klein et al.|[2012|l. The central position of 

oil ,,,00 o,_— x _— 


all observations was q'[j2ooo]= 12"01'"36?3, h[j2ooo]=-65‘’08'53'.'0. 
We checked the focus at the beginning of each observing ses¬ 
sion, after sunrise and sunset on Saturn. We checked line and 


' This publication is partly based on data acquired with the Ata¬ 
cama Pathfinder Experiment (APEX). APEX is a collaboration be¬ 
tween the Max-Planck-Institut fur Radioastronomie (MPIfR), the Eu¬ 
ropean Southern Observatory, and the Onsala Space Observatory. 
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Fig. 1. The BHR71 bipolar outflow in its entirety, in Left: CO (6-5) (colours, associated with the wedge, in main beam temperature units, K km s“*) 
and (3-2) (white contours, from 30% to 100% of the signal in steps of 10%), both observed by the APEX telescope and integrated between -50 
and +50 km s“*, with the resolutions indicated by green circles in the lower right comer; Right: the 8 gim emission detected by the SpitzerflKAC 
receiver (colours, N09), with the H 2 0-0 S(5) emission as observed by the SpitzerflRS receiver in the inner parts of the outflow (white contours, 
from Gil), and the SiO (5^) emission (red and black contours, Gll) in the upper lobe (the green circles in the lower right comer show the 
respective resolutions of CO (3-2), (11-10), and SiO (5^)). On both maps, the grey inset is the field shown in figure]^ the SiO peak and knot 
positions are indicated in blue, the knot 5 and HH320 region are in pink, and the IRS 1&2, knots 6, 8, and HH321 region are indicated by grey 
dots. 


continuum pointing locally on IRC+10216, 07454-7112, or t] 
Car. The pointing accuracy was better than ~5" rms, regard¬ 
less of which receiver we used. Table [T] contains the main 
characteristics of the observed lines and corresponding ob¬ 
serving set-ups. The observations were performed in position¬ 


switching/on-the-fly mode using the APECS software ( |Muders| 
|et al.||2UD6l l. We reduced the data with the CLASS software 
(see http;//www.iram.fr/IRAMFR/GILDAS). This reduction in¬ 
cluded baseline subtraction, spatial, and spectral regridding. For 
all observations, the maximum number of channels available in 
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Fig. 2. Overlay of the velocity-integrated CO (6-5) (colour background) 
with the CO (3-2) (white contours) emission observed with the APEX 
telescope in the inner upper lobe of the BHR71 outflow. For both lines, 
the intensity was integrated between -50 and 50 km s“*. The wedge 
unit is K km s“* in main beam temperature. The CO (3-2) contours are 
from 50 to 100% of the maximum, in steps of 10%. The half-maximum 
contours of the CO (3-2) and (6-5) maps are indicated in red and black, 
respectively. The dark blue circles indicate the positions and beam size 
of the SOFIA/GREAT observations for the CO (11-10) transition. The 
APEX and SOFIA beam sizes of our CO (6-5), (16-15) and (11-10) 
observations are also provided (upper left comer light green circles, see 
also table[^. The pink hexagons mark the position of the so-called knot 
5 and HH320AB object. 



a[J2000] 


the backend was used (8192). We obtained maps for all consid¬ 
ered transitions, covering the held of the whole outflow shown 
in flgure[T] 

Figure 1^ shows the velocity-integrated CO (6-5) broad-line 
emission (colours, resolution 9") in the upper inner part of 
BHR71, overlaid with the CO (3-2) (white contours, resolution 
18")- The higher angular resolution of the CO (6-5) line emis¬ 
sion shows that it traces the walls of the cavity of the outflow 
associated with IRSl (see Figure [^l, with a local maximum of 
emission also corresponding to the outflow driven by IRS2 (the 
HH320 AB position in the map). The outflow is seen close to the 
plane of the sky. The two positions of interest identified by G11 
are marked by a blue hexagon and circle corresponding to the 
beam of our CO (11-10) observations with GREAT (which will 
be the focus of our analysis, see sections [23] and[^. The position 
of the SiO peak, detected through 28"-resolution observations of 
the SiO (5-4) is localised between two emission peaks in CO (6- 
5). The most prominent of these CO peaks is the southern one, 
which corresponds to the so-called SiO knot; it coincides with 
an H 2 emission peak (also see figure]^ which overlays the same 
CO (6-5) held with the H 2 0-0 S(5) data of N09 and Giall). 
The half-maximum emission contours are also given in this map 
in thick black and red contours for the (6-5) and (3-2) lines, 
respectively, at their nominal resolutions. 
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2.3. SOFIA-GREAT observations of CO 


The observations towa rds BHR71 were con ducted with the 
GREAlj^ spectrometer ( Heyminck et al.|[2012 i during SOFIA’S 
Cycle 1 ‘southern deployment’ flight on July 23rd, 2013. We 
observed two positions, towards the northern apex of the inner 
outflow structure; the SiO knot and peak, as defined in G11 (fig¬ 
ures [T] and [^. We tuned the receiver to the CO (11-10) and 
(16-15) lines frequency 1267.014 GHz LSB and 1841.346 GHz 
USB. We connected the receiver to a digital FFT spectrometer 
( Klein et dr||2012 1 providing a bandwidth of 1.5 GHz with re¬ 
spective spectral resolutions of 0.018 and 0.012 km s '. The ob¬ 
servations were performed in double beam-switching mode, with 
an amplitude of 80" (or a throw of 160") at the position angle of 
135'’ (NE-SW) and a phase time of 0.5 sec. The nominal focus 
position was updated regularly against temperature drifts of the 
telescope structure. The pointing was established with the optical 
guide cameras to an accuracy of ~5". The line and observation 
parameters are listed in table The integration time was 13 min 
ON source for each line, for respective final r.m.s of ~0.70 and 
0.75 K. The data were calibrated with the KOSMA/GREAT cal¬ 
ibrator ( |Guan et al.|2012| l, removing residual telluric lines, and 
subsequently processed with the CLASS softwar^ 


3. Results: the CO data 

3.1. Spectra 

Figure [^presents all '^CO spectra obtained in the two targeted 
positions in the course of our APEX (Tup = 3,4,6,7) and SOFIA 
(Tup = 11, 16) observations, after convolution to the same angu¬ 
lar resolution, that of the CO (11-10) transition (see table[T]for 
observing parameters related to each of these lines). The only ex¬ 
ception to this convolution is the CO (16-15) line, for which only 
a single-point observations is available. As the line was observed 
at an angular resolution of roughly 1673, we could in principle 
use the integrated intensity extracted from this line as an upper 
limit to that convolved to the 24" resolution. In fact, we cor¬ 
rected this value for the beam dilution effect in our analysis, see 
discussion in section lSTI on how we made use of this line. On the 
figure, we split the spectra in two groups for visibility purposes 
(Tup= 3,4,6 in the upper panels, and 6,1, 11, 16 in the lower 
panels). 

In both positions, all lines exhibit a profile typical of the pres¬ 
ence of a pure shock, with wings extending towards red-shifted 
velocities, up to 30^0 km s“*; up to Tup = 6, self-absorption is 
also detected at the velocity of the cloud (around -4.5 km s '). 
The CO (3-2) and (4-3) profiles coincide very well, suggesting 
that these lines are optically thick. This is confirmed for the CO 
(3-2) line, for which we observed the *^CO isotopologue on both 
positions (see section [3T2| table[T]and figurej^for the spectra ob¬ 
tained on both positions). In both positions, the CO (6-5) profiles 
start to differentiate from the lower-lying profiles. This departure 
is confirmed in the higher-lying transitions, and reveals that the 
lines from Tup = 6 are probably excited in different layers of the 
shocked region than their lower-lying counterparts. Globally, it 
can be seen that the excitation conditions vary with the posi¬ 
tion: the Tup = 11 and 16 profiles differ, and the CO (16-15) is 

^ GREAT is a development by the MPI fur Radioastronomie and the 
KOSMA/Universitat zu Kbln (Kdlner Observatorium fur SubMillime- 
ter Astronomie), in cooperation with the Max-Planck-Institut fur Son- 
nensystemforschung and the Deutsches Zentrum fur Luft- und Raum- 
fahrt Institut fur Planetenforschung. 

3 http://www.iram.fr/IRAMFR/GILDAS 
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Table 1. Observed lines and corresponding telescope parameters for the APEX and SOFIA observations of BHR7I. 


species 

bJ 

n 

o 

n 

o 

f3 

n 

o 


n 

o 


‘"CO 

line 

(3-2) 

(^3) 

(6-5) 

(7-6) 

(11-10) 

(16-15) 

(3-2) 

telescope 

APEX 

APEX 

APEX 

APEX 

SOEIA 

SOEIA 

APEX 

y (GHz) 

345.796 

461.041 

691.473 

806.652 

1267.014 

1841.346 

330.588 

FWHM (") 

18.1 

13.5 

9.0 

7.7 

24 

16.3 

18.9 

sampling (") 

6 

6 

4 

4 

pointed 

pointed 

6 

receiver 

FLASH345 

FLASH460 

CHAMP+ 

CHAMP+ 

GREAT 

GREAT 

FLASH345 

observing days 

2013-06-03 

2013-06-03 

2013-06-04 

2013-06-04 

2013-07-23 

2013-07-23 

2013-06-03 


2013-06-04 

2013-06-04 

2014-06-28 

2014-06-28 

2013-07-23 

2013-07-23 

2013-06-04 

CefF 

0.95 

0.95 

0.95 

0.95 

0.97 

0.97 

0.95 

Beff 

0.69 

0.60 

0.42 

0.35 

0.67 

0.67 

0.69 

Tsys (K) 

169-193 

408-610 

3650-10000 

12000-34000 

3274-3307 

3160-3207 

182-224 

An (km s“') 

0.331 

0.496 

0.635 

0.544 

0.018 

0.012 

0.346 

reference offset (") 

(-120,260) 

(-120,260) 

(-120,260) 

(-120,260) 

beams witch 

beamswitch 

(-120,260) 



Fig. 3. CO transitions in the SiO peak (left panels) and knot (right panels) positions indicated in figure]^ APEX (3-2), black line ; (4-3), red line ; 
(6-5), green line and histograms (upper panels); (7-6), dark blue lines ; SOFIA (II-IO), pink lines and (16-15), grey lines (lower panels, overlaid 
on the green histograms of the 6-5 transition). The last two were multiplied by five for comparison purpo ses. Respective spect ral resolutions are 


0.33, 0.50, 0.64, 0.54, 0.90 and 0.62 km s *. The vertical dotted line marks the cloud velocity, -4.5 km s * (Bourke et al. 


1995b I. 


not even detected in the SiO peak position. For both positions, 
the excitation conditions also vary with the velocity: close to 
the systemic velocity, the lower-7 line emission greatly domi¬ 


nates, an effect previously reported in LI 157 B1 by Lefloch et al. 
( |MT2t . Because of the associated weak or absent detection in the 
SOFIA lines, and because the two selected positions are not in¬ 
dependent, we decided to exclude the SiO peak from the shock 
analysis presented in section]^ (also see an additional argument 
in section|4. Ik 


3.2. CO (3-2) line opacities 

Figure shows our spectra of the '^CO and '^CO lines in the 
SiO peak and knot positions at the same spectral and spatial 
resolutions. These spectra allow us to plot the ratio of the line 
temperature. For each velocity channel, the line temperature ra¬ 


tio '^CO / ‘^CO is hence also shown (right ordinates) with the 
following colour-code: red dots show the ratios for the velocity 
channels where the '^CO is detected at more than 3cr, and orange 
dots where the detection was obtained with a confidence 
level between 2 and 3cr. The grey dots are for all the other chan¬ 
nels. Horizontal dashed lines show the reference values of 20 
and 40 for these ratios. The orange and red dots all lie below 30. 
Assuming an identical excitation temperature for the '^CO and 
'^CO lines, and a typical interstellar abundance ratio of 50-60 
(e.g. Langer & Penzias||1993[) , these line ratios yield minimum 
opacity values of 3 for the ‘^CO lines. We therefore conclude 
that the emission is optically thick at least in the low-velocity 
regime of the spectral wings. The large optical thickness is con¬ 
sistent with the constancy of the line integrated intensities (in 
the non-absorbed components) from CO (3-2) and (4-3), when 
convolved to the same resolution(s). As our observations yield a 
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Fig. 4. Comparison of the '^CO (black line) and '^CO (pink line) (3-2) 
emission profiles obtained in the SiO peak (upper panel) and knot (lower 
panel) positions. The ratio of main-beam temperatures is also shown in 
the form of red dots (for the channels where the '^CO detection is over 
3cr), orange dots (for the channels where the '^CO detection is between 
2 and 3cr), and grey dots (for the remaining channels). Associated with 
this distribution are the grey horizontal dotted lines, which show values 
of 20 and 40 for this ratio. 


minimum optical thickness of 3 up in the wings, we decided to 
exclude both the CO (3-2) and (4—3) lines from our analysis. 


3.3. Dynamical age of the outflow 

An important parameter for the modelling of young outflows is 
their age. As we have mapped the entire outflow, we are able 
to give an upper limit of its age, simply obtained by dividing 
the distance between the furthest point from the driving proto¬ 
star and the protostar by the associated linewidth. The positions 
we consider belong to the northern lobe of the outflow powered 
by the IRS 1 protostar. As the furthest point with significant CO 
emission belonging to this outflow, we adopted the furthest point 
on the 10% CO (3-2) emission contour or the furthest point with 
3cr emission away from the driving protostar, and found that 
these two approaches were equivalent. Indeed, the offset from 
both points relative to the IRSl position is about (-72", 280"), 
corresponding to a distance of ~289". At a distance of 200 pc, 
this translates in ~5.8xl0'* AU. In this position, the full-width at 
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Fig. 5. Overlay of the map of CO (6-5) emission observed by the APEX 
telescope (colour background) with the H 2 0-0 S(5) emission (white 
contours), observed with the Spitzer telescope. The wedge unit is K km 
s“* (main beam temperature) and refers to the CO observations. The 
H 2 0-0 S(5) contours are from 10% to 100%, in steps of 10%. The light 
blue contour defines the half-maximum contour of this transition. Like 
in figure the blue circles and dots indicate the positions and beam 
sizes of the SOFIA/GREAT observations. The beam and pixel sizes of 
the CO (6-5) (11-10), (16-15), and H 2 0-0 S(5) observations are the 
green circles in the upper left corner. The black contour delineates the 
half-maximum contour of the Fl 2 0-0 S(2) transition. The field is the 
same as in figure]^ and the knot 5 and FIFI320 region are in pink. The 
field covered by Sp/tzer/TRS to observe the H 2 emission is indicated in 
dashed grey line. It excludes the SiO peak position. 



zero intensity of both the CO (3-2) and (6-5) (unshown) lines 
is about 15-20 km s“\ which converts to a dynamical age of 
(1.4-1.8)xl0"^ years. Of course, this is an upper limit of the real 
age of the outflow, as it is likely that the apex of the outflow was 
associated with greater velocities in the past. A similar measure¬ 
ment for the SiO knot position (i.e. dividing its distance to the 
protostar by the full-width at zero intensity) yields a dynamical 
age of ~1800 years. 


4. Results: existing data 

4.1. H 2 observations 

We used the Spitzer/IRS observations of the H 2 pure rota¬ 
tional transitions, 0-0 S(0) up to S(7), reported and analysed 
in N09, Giall, and Gil. The reduced data kindly communi¬ 
cated to us by David Neufeld contain rotational transition maps, 
with 3'.'6 angular resolution, centred on a[j 2 ooo] = 12^'01'"36?31, 
<5[J2000]=-65°08'53'.'02. Figure|^shows an overlay of our APEX 
CO (6-5) map with the H 2 0-0 S(5) region observed by Spitzer, 
both at their nominal resolution. The figure shows coinciding 
maxima between the two datasets in the selected position, and a 
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slightly different emission distribution. This might be the effect 
of the better spatial resolution of the H 2 data, which reveal more 
peaks than in CO. This overlay also shows the similar morphol¬ 
ogy of the emission of the S(2) (half-maximum contour in black) 
and S(5) (half-maximum contour in light blue) transitions on the 
region we chose to analyse. Unfortunately, the hgure also shows 
the coverage of the H 2 emission observations, and that the SiO 
peak position was not covered by the H 2 observations, which is 
another reason to exclude it from our shock analysis presented 
in section |5] 

As part of our modelling (see Sect.|^, we used the excitation 
diagram derived for the selected emission region around the SiO 
knot position. The H 2 excitation diagram displays \n{N„j/gj) as 
a function of E„jlk^, where N„j (cm“^) is the column density of 
the rovibrational level (v, 7), E^jlks is its excitation energy (in 
K), and gj - (2] + l)(2/-i-1) its statistical weight (with 7=1 and 
7 = 0 in the respective cases of ortho- and para-H 2 ). If the gas 
is thermalised at a single temperature, all points in the diagram 
thus fall on a straight line. 

The initial resolution of the maps is 3'.'6, but we wanted to 
operate at the same resolution as for CO. Unfortunately, the SiO 
knot position lies at the edge of the field covered in the H 2 map, 
as can be seen in figure]^ We consequently used three different 
methods for extracting the H 2 line fluxes. First, we extracted the 
flux from the initial map at its nominal 3'.'6 resolution; second, 
we used the first method, but applied to a map that was convolved 
to the 24" resolution; and third, we associated a filling factor of 
0.2 to the fluxes obtained through the second method. 

Because of the location of the SiO knot and the rather large 
beam of our analysis, the first two methods provide lower limits 
to the H 2 line fluxes. On the contrary, with method (3) we volun¬ 
tarily extracted upper limits to these fluxes. We then used an av¬ 
erage value between the values inferred from the three methods, 
and computed rather large errorbars based on the combination 
of all three methods. The values were corrected from interstellar 
extinction following the treatment already applied in Gil. The 
result can be seen in the excitation diagram shown in the lower 
panel of figure]^ where points corresponding to methods (1), (2), 
(3) are shown in black empty diamonds, upward and downward 
triangle, and the resulting average points are the black squares, 
respectively. The overall average H 2 values we used to build our 
figure are given in table|B.l| 


4.2. SiO 

As part of our study, we also decided to re-analyse the SiO ob¬ 
servations already presented in Gil. The emission from three 
lines was mapped in the northern lobe of the outflow; SiO (5-4), 
(6-5), and (8-7), and the corresponding spectra were extracted 
in the SiO peak and knot positions. As maps were obtained, 
all spectra were convolved to the resolution of the SiO (5^) 
line, about 28". In the present study, we hence corrected the SiO 
dataset for the slight difference in resolution between this value 
and that of the CO (11-10) observations, 24", by simply mul¬ 
tiplying the Gil integrated intensities by (28/24)^. Also, given 
the rather large beam size, we generated an integrated intensity 
diagram (displaying J TmecIc against J^p for each transition) for 
three different filling factor values, 1, 0.5, and 0.2. The result can 
be seen in the form of respective black, dark grey and light grey 
squares in figure [7] T he SiO values we used to build our figure 
are listed in table lB.ll 


5. Discussion 

In the following we seek to constrain the physical conditions in 
the shocked regions. Given the large number of observational 
constraints, covering several species and several associated tran¬ 
sitions, the method-of-choice for obtaining these physical con¬ 
ditions is a comparison to sophisticated shock models. The re¬ 
sults of the shock modelling forms the foundation of the discus¬ 
sion; what characterises the chemistry of the outflow, in particu¬ 
lar with respect to SiO, and what characterises the energetics. 


5.1. Constraining shock modeis 


Shock models are used to constrain the physical conditions in the 
outflow shocks through comparison with H 2 and CO, following 
the methods already used in [Gusdorf et al. ( |2008a| |2012| l and 
Anderl et al. (2014|l. We generated CO flux diagrams and H 2 ex¬ 


citation diagrams from the observations of the SiO knot position, 
and compared them with results from a grid of models computed 
with the ‘Paris-Durham’ (e.g. [Flower & Pineau des Forets|2003] l 
1-D shock model. The H 2 diagram has already been introduced 
in section [4.1 j and is shown in figure Concerning the observ¬ 
able associated with CO, we chose to display the data in the form 
of a flux diagram (line flux vs. J^p), because of the increasing 
number of extragalactic studies that use this form. In the present 
case, we benefit from the fact that our observations were pointed 
towards a pure shock position. To generate a CO flux diagram, 
we consequently integrated the signal obtained over the whole 
line profile at the angular resolution of 24", associated with a fill¬ 
ing factor of 1 (see sectionj^for an explanation on filling factors 
for various CO transitions). Because of their opacity, and their 
likeliness to be contaminated by ambient gas, we excluded the 
(3-2) lines from our analysis (see section [T2| i. Given the opacity 
values inferred from this line, we also excluded the (4-3) line 
from our fitting procedure. Overall, the two observational tools 
we used can be seen in figure a flux diagram for CO (upper 
panel) and an excitation diagram for H 2 (lower panel) on which 
the observational points for the SiO knot are always shown in 
black squares. The CO values we used to build our figures are 
listed in table EH The observational flux diagram for the SiO 
peak position can be found in Section |C| (Figure [C.lj l. 

The grid of shock models is that of G11, also used in other 
publications already mentioned ( [Gusdorf et al.||2()12[ [Anderl[ 
[et al.[2014] ). To summarise, we used a grid covering the following 
input parameters; pre-shock density tin from 10^ to 10^ cm“^, 
magnetic field parameter from b - 0.3 to 2 (defined by B(p G) = 
b X [mh (cm“^)]'^^, where B is the intensity of the magnetic field 
transverse to the shock direction of propagation), and shock ve¬ 
locity Vs from 15 to 35 km s '. Our grid contains both stationary 
C- and J-type models, and also non-stationary, so-called CJ-type 
models ( [Lesaffre et al.|[2004a|b| l. A more complete description 
of the parameters space coverage {tin, b, and Vs) can be found in 
Table The parameter coverage is not complete; not all veloc¬ 
ities are present in our grid for all values of the magnetic-field 
parameter, b. In fact, the velocity of C-type shocks must remain 
below a critical value that depends mainly on the pre-shock den¬ 
sity and magnetic-field parameter (Le Bourlot et al.|2()02 Flower 


[& Pineau des Forits][2003[ ), which explains the decrease of the 
maximum shock velocity with the pre-shock density in our grid. 
Additionally, the time necessary to reach a steadystate depends 
on the pre-shock density, shock velocity, and magnetic field val¬ 
ues. Following the method presented in G08b (Section 4.1), the 
set of C-type shock models enabled us to restrict the range of the 
search in the parameter space for the CJ-type shock models. We 
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Table 2. Shock model parameters. Grid intervals Ax are listed as minimum and maximum increments found in the grid. 


Shock type 

Number of models 

V [km s ^ ] 

Av [km s 

b 

Ab 

wh [cm 

Age [yr] AAge [yr] 

C-type" 

108 

20-55 

2-5 

0.45 - 2 

0.15-0.25 

10^ 10^, 10\ 10'’ 

_ _ 

C-type* 

9 

20-40 

10 

1-3 

0.5 

10^ 

- - 

J-type 

21 

10-50 

2-15 

0.1 

- 

10'', 10^ 

- 

CJ-type 

1035 

10-50 

2-5 

0 . 3-2 

0.15-0.25 

10^ 10^ 10^ 

4- 15000 1 - 1285 


Gusdorf et al. 

2008b 

2011 

2012 1 '*' 

Anderl et al. 

2014 


then computed a grid of non-stationary shock models around a 
first estimate of the shock age, making sure that the range of ages 
was sufficient to include any model likely to fit the H 2 observa¬ 
tional data. The final considered ages range from a few tens to 
around fifteen thousand years. Our grid also includes a few sta¬ 
tionary models including the effects of grain-grain interactions. 


B1 (e.g. Gusdorf et al. 2008b Flower & Pineau des Forets 2012| l. 
It is also the value associated with our ten best fits to the dataset; 
it is relatively well constrained, as in particular the general shape 
of the CO diagram is very sensitive to the density, i.e. to the pre¬ 
shock density parameter. In our ranking of models by decreasing 
value, the first model with a different pre-shock density has a 


as presented in Guillet et al. 

(201 ip, Anderl et al. (2013|>, and ^ value that is more than hve times the lowest one. Generally 

already used injAnderl et al.|( 

2014|l,|Leurini et al.|(|2014a|l. speaking, it is a very standard value when modelling shocks from 


We compared models and observations based on the three 
higher-7 CO lines, from CO (6-5) to (11-10), and on all the 
H 2 pure rotational lines. The CO (3-2) and (4-3) lines were ex¬ 
cluded from our procedure, as stated above. Our initial CO (16- 
15) observation is in principle an upper limit, since it was as¬ 
sociated with a slightly smaller beam than the other lines. Nev¬ 
ertheless we assumed that the CO (16-15) emission is as ex¬ 
tended as the CO (6-5) emission, and we estimated the resulting 
flux over a 24" beam by multiplying the value observed over a 
16'.'3 beam by a factor (16.3/24)^. We then treated this corrected 
value as any other CO observation. Similar to what we did in 
previous studies, we used a routine to compare these obser¬ 
vations to the models. The best results can be seen in figure 
On each of these panels, we show the points of our best fit in 
red circles, plus three other satisfying models in smaller, purple, 
green, and blue circles. Our best-fit model (in red points) is non- 
stationary, with the following characteristics: riu - 10^ cm“^, 
b = 1.5, Us = 22 km s“', and an age of 3800 years. The;^f^ value 
of our 10 and 20 best-fit models are within a factor 1.4 and 2.1 to 
the best (lowest) value. Broadly speaking, we found the H 2 lines 
in particular very difficult to fit: tfie exact curvature of tfie exci¬ 
tation diagram is only approacfied by our best model. We tried 
in vain to improve the quality of the fit by changing the ortho-to- 
para ratio of H 2 in our calculations (switching its value from 3 
to 1 and 2). Eventually, we can not exclude that a better fit could 
be found via a better gridding of our parameters space. However, 
it turned out that this best-fit model also fitted the CO (3-2) and 
(4-3) lines quite in a very satisfying way. These values can be 
compared to what was constrained by G11: hh = 10"^ cm“^, b - 
1 - 2, Us = 25 -30 km s“', and an age of 300-800 years. A few 
comments can be made on these shock parameters. 

Shock type. The shock type we constrained is the same as 
found by Gil. Indeed, it is a very general result that in low- 
mass star forming environments, H 2 excitation diagrams can be 
fitted by these kinds of models only (e.g. |Giannini et al.||2006 


[Gusdorf et al.|2008bl [Flower & Pineau des Forets|2013| l. All of 
our ten best fits are CJ-type models. Even more, the 7- and C- 
type models are all at the end of our fist of best-fit models: the C- 
type models generally do not fit the pure rotational H 2 excitation 
diagram, while the 7-type models do not generate enough CO 
emission to match the observations. 

Pre-shock density. Eirst, the pre-shock density value is the 
same as that found in previous analyses of BHR71 (Giannini 
et al.|2004 for the analysis of the HH320 region, Gil for the SiO 


knot position), and in the resembling pure shock position LI 157- 


was also found in shocks associated with SNRs that are interact¬ 
ing with the interstellar medium ( jCesarsky et al.|1999| |Gusdorf| 
jet al.|201^|Anderl et al.|20l4] l. This corresponds to a post-shock 
density of ~10^ cm a value that is one or two orders of mag¬ 
nitude smaller than those found by Gial 1 in the HH320 and 321 
regions of the same outflow. This could be due to the fact that 
their study only focussed on the brightest pixel associated with 
these regions, whereas we consider a rather large beam size. On 
the other hand, our value is just above the H 2 density averaged 
over the whole inner outflow found by N09 (10^'* cm"^). This 
can be explained by the fact we are focussing on a bright H 2 
spot, and also that the post-shock value we are providing here is 
that of a stationary post-shock. In the (realistic) case where the 
shock has a finite spatial size, tfie post-sfiock gas will expand and 
come into pressure equilibrium with its surroundings, resulting 
in a globally lower value than ours. 

Magnetic field strength. The strength of the transverse mag¬ 
netic field is 150 pG in tfie pre-shock region, that turns into 
1.62 mG in the post-shock, given the compression factor and 
the fact that the magnetic field is frozen in the neutral fluid. It is 
comparable to what was constrained in Gil {b - 1-2), and also 
to what was found in the studies of HH320 in BHR71, LI 157- 
Bl, or HH54, or in the aforementioned SNR studies. This value 
is also consistent with the analysis of Zeeman measurements in 
molecular clouds where magnetic and kinetic energy densities 
are in approximate equipartition ( |Crutcher]| 1999) 1. Our ten best 
fits all predict a value of the b parameter between 1.0 and 2.0. 
Eigure shows the best model with a b value outside of this 
range, in purple points. This model is associated to a x^ value 
that is 1.8 times the lowest. This model is also a bit older than 
our ten best-fit ones (see below). It does not fit the lower-7 data 
as satisfyingly as our best-fit model (red points in Eigure]^. 

Shock velocity. Tfie shock velocity values of our ten best- 
fit models all are between 20 and 25 km s '. In our ranking of 
models by decreasing x^ value, the first model witfi a velocity 
out of this range has ?lx^ value that is about three times the low¬ 
est one. More specifically, tfie sfiock velocity tfiat we constrain is 
slower tfian, e.g. tfie full-widtfi at zero intensity of tfie CO lines. 
The modelled shock velocity is the difference between the pre¬ 
shock velocity and the impact velocity. Consequently, if the pre¬ 
shock material is already accelerated, the actual shock velocity 
is expected to be correspondingly lower. Preliminary accelera¬ 
tion of the pre-shock is possible in the SiO knot position, as it 
has already been encompassed by the shock that is now propa- 
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Fig. 6 . Our model-observations comparisons. Upper panel: CO flux di¬ 
agram over a beam of 24"; the observations in the SiO knot position 
are the black squares, and the model results are the coloured circles (see 
text). The CO (16-15) observational point is corrected for beam-size 
effect. The uncorrected point is the upper limit (see text), indicated by 
the grey arrow in this panel. Lower panel: H 2 excitation diagram for 
the SiO knot position, extracted following the different procedures de¬ 
scribed in the text (empty symbols), global average in black squares, 
and model results in coloured circles (see text, with the same code as in 
the upper panel). 


gating at the northern tip of the outflow (a bit more than twice 
as far from the driving protostars). Alternatively or simultane¬ 
ously, this velocity discrepancy could be the sign of a limitation 
to our models. Indeed, we are modelling a multi-dimensional 
complex shock region, which is seen edge-on, by means of a 
one-dimension model, which is considered to be face-on. Addi¬ 
tionally, it is likely that the rather large beam of our observations 
contains not one shock structure, but a collection of them. This 
can be seen in the spatial sub-structure of the H 2 emission re¬ 
vealed in figure where multiple peaks are detected. This can 
also be seen in the CO (figuresandor SiO (Gil) line pro¬ 
files, where a ‘plateau’ or a ‘bump’ can be seen between 20 and 


40 km s ', very much resembling the bullets revealed in the CO 
(in e.g. Cep E, Gomez-Ruiz et al.|2012|l or water line profiles (in 
e.g. L1448-mm, [Kristensen et al.|2011| l. Another expression of 
the intrinsic multi-dimensional nature of the observed shocks is 
the impossibility of fitting the different molecules with the same 
filling factor values. For instance, H 2 is only associated with the 
hottest regions in our beam, as a 500 K temperature is necessary 
to populate its levels, whereas CO is more easily excited (for 
7up = 6, Eup = 116 K). In particular, unlike H 2 , CO is probably 
emitting in the post-shock expansion region mentioned in the 
‘pre-shock density’ paragraph above, where the primary shock 
structure progressively mixes with the ambient medium in all 
spatial dimensions. This effect cannot be accounted for by a 1-D 
shock model. 

Shock age. Finally the shock age is rather large; it is larger 
than what was previously constrained for this region of BHR71 
by Gll (300-800 years), which is due to the high amounts of 
CO that we detected, that are not compatible with younger CJ- 
type shocks. However, again the age predicted by our ten best- 
fit models consistently lies between 3500 and 5500 yr, with the 
model in purple in the Figure l^being slightly older. In our rank¬ 
ing of models by decreasing ;^"^value, the first model with an age 
younger than 3000 years has a value that is about twice the 
lowest. The age that we find is between the dynamical age of the 
SiO knot and that of the global northern lobe (see Section [33] l, 
which hints again at the fact that we are probably catching sev¬ 
eral shock episodes in our beam. 

Limitations. The question of the largest source of uncertainty 
in these results is difficult to assess. The observational problems 
(the different beam sizes and the fact that the SiO knot lies close 
to the edge of the H 2 observations) make for an important lim¬ 
itation, but we believe we have taken it into account by using 
conservative errorbars. More problematic is the complex nature 
of shocks structures. Lacking high angular resolution, we could 
only assume that we catch a single shock structure in our beam. 
Our results seem to indicate this is not the case, given the con¬ 
strained shock velocity and ages, and the different filling factors 
we had to adopt for each molecule. A shock similar to those 
we have in our grid of models is probably propagating in one 
beam of 24", but it is very likely to be associated with less vi¬ 
olent shocks corresponding to the mixing of its warm, dense, 
and accelerated post-shock with the surrounding ISM. To sum¬ 
marise, we could only be aware and accept the fact that we are 
averaging processes out by working over a beam size typical of 
single-dish observations. From the modelling point of view, so¬ 
lutions do exist to more realistically account for the observations 
of complex geometrical shock structures. The first consists of 
adopting a probability density function (PDF) of shocks, which 
is to consider that the observed shock characteristics follow a 
statistical distribution. This distribution can unfortunately only 
be guessed, and hence generates additional free parameters (see 
[Lesaffre et al.|20T3] for applications of this method). The sec¬ 
ond consists of stitching 1-D shock layers (similar to the one we 
consider here) to either a curve (e.g. Kristensen et al.|2008) l or a 
bow-shock structure (e.g. |Gustafsson et al.|2010| l, hence simulat- 
ing ‘pseudo-’ 2- or 3-D shock propagation. This approach also 
yields free parameters, such as the inclination of the magnetic 
field with respect to the shock structure. None of these methods 
would be ideal in the present case, given the lack of knowledge 
of the small-scale shock structure within the 24" beam (no in¬ 
dication on the collection of shocks is available at the moment), 
and the relatively limited dataset (if more free parameters are to 
be considered, more constraints are needed to lift the degener¬ 
acy in the results). Their application will become relevant when 
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Fig. 7. The integrated intensity diagram generated from the G1 1 obser¬ 
vations of three SiO lines (Tup = 5, 6, 8), corrected for beam size and for 
filling factor effect: filling factor of 1, 0.5, and 0.2 in black, dark grey, 
and light grey squares. The result of our best-fit model for the H 2 and 
CO lines (uh = 10"* cm“^, b = 1.5, = 22 km s“*, and an age of 3800 

years) is shown with 1, 2, or 4% of the pre-shock Si placed in the grain 
mantles in red empty, dotted, or filled circles. 


the SiO knot is observed at higher angular resolution (in CO and 
H 2 ), and if possible with a maximum number of observations in 
key species (such as H 2 , CO, and SiO, but also H 2 O, OH and O i 
for instance). We note that at this point, a more tightly sampled 
shock model will prove necessary. 

In the following we will present results based on our best 
fit only. All of our ten best-fit models were found to be non¬ 
stationary (C7-type), with hh = 10"^ cm“^, b — 1 - 2, Uj = 20 - 
25 km s“', and an age of 3500^500 years. 


5.2. Employing shock models: SiO chemistry 


After constraining shock models based on CO and H 2 obser¬ 
vations, it is now possible to fine-tune the SiO chemistry. SiO 
has long been interpreted as a tracer of C -type shocks with a 
velocity greater than 20-25 km s“' (e.g. 


Caselli et al. 


1997 


ISchiUce et al.|1997{|Gusdorf et al.|2008a| l. Indeed, in these types 
of shocks, the drift velocity between the charged grains and the 
most abundant, neutral molecular species is sufficient to gener¬ 
ate the sputtering of the core of the grains, where all the silicon¬ 


bearing material was considered to be locked. However Gusdorf 
|et al.| ( (2()08b| l have demonstrated: 1) that CJ-type shock mod¬ 
els were the only ones to fit the H 2 emission in L1157-B1, and 
2 ) that this process was not efficient enough to generate levels 
of SiO emission comparable to the observations in this kind of 
(CJ-type) shock models. In order to be able to self-consistently 
account for both SiO and H 2 emission, one solution was consid¬ 
ered, consisting of transferring a fraction of SiO from the core to 
the mantle of the grains in the pre-shock phase. The maximum 
fraction of SiO thus transferred to the mantles was set to 10%. 
With mantle sputtering being easier than core sputtering, these 
authors were then able to fit the H 2 and SiO emission by means 
of a single, non-stationary shock model. The presence of silicon¬ 
bearing material in the grain mantles has also been assumed in 
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ICoutens et al.| ( |2OT3] l to explain the narrow emission component 
detected in the SiO (2-1) line profile at the systemic velocity 
around NGC1333 IRAS 4A. 

Interestingly, the transfer of Si from the core to the mantles 
in the form of SiO was also considered by |Anderl et al.| ( [2()13| l in 
the context of denser shock models. Indeed, for pre-shock densi¬ 
ties above 10 ^ cm“^, the effect of grain-grain interactions cannot 
be neglected in shock models, as shown by |Guillet et al.| ( |201 l| l. 
Taking these interactions into account results in the significant 
production of small, charged grains that couple very well with 
the neutral fluid, in effect reducing the width of the shock layer, 
and increasing its maximum temperature. In the narrow shock 
layers thus produced, Anderl et al. ( |2013[ ) have shown that the 
grain core sputtering is not efficient enough to produce levels of 
SiO emission comparable to the observations, hence the recourse 
to the transfer of a fraction of Si towards the mantles in the pre¬ 
shock phase. [Leurini et aT] ( |2014a] l have validated this approach 
by successfully confronting these models with observations. 

Consequently, we ran our best-fit model for the H 2 and CO 
data(nH = 10"^cm~^,h= 1.5, Us = 22 km s“', and an age of 3800 
years), including 0 to 10 % of the pre-shock silicon-bearing ma¬ 
terial in the form of SiO in the grain mantles. We then generated 
the corresponding integrated intensity diagrams, which we com¬ 
pared with the observations. The result can be seen in figure |7] 
For a filling factor value conservatively varied between 0.2 and 
1, we show that no more than 4% of SiO needs to be placed in 
the grain mantles to reproduce the observations. The presence of 
silicon-bearing material on the grain mantles could be explained 
by the fact that the SiO knot position has already been processed 
in the past by the shocks that are now located at the apex of the 
northern outflow lobe, further north, as can be seen in figure 
Under this assumption, we have demonstrated that it is possi¬ 
ble to self-consistently fit H 2 , and velocity-resolved CO and SiO 
observations in the ‘SiO knot’, pure shock position of BHR71. 


5.3. Employing shock models: energetics 

In this section, we discuss the energetics associated with the 
shocks along two axes: the CO excitation generated from our 
best-fit model, and the derivation of the outflow parameters. 

CO excitation. Indeed, the excitation of CO has been inten¬ 
sively used in the literature to demonstrate the presence of vari¬ 
ous processes at work in the regions of star formation. For exam¬ 
ple, [yan^Smpenjet^l 2010^ and ||Visser et al.| ( |2012] l obtained 
CO flux diagrams (similar to our figure with PACS, centred 
on low-mass protostars at the origin of various outflows (HH46, 
NGC1333 IRAS 2A, and DK Cha). Based on these diagrams, 
they evidenced the existence of three physical components cor¬ 
responding to three distinct processes contributing to the excita¬ 
tion of CO; a passively heated envelope, UV-heated outflow cav¬ 
ity walls, and small-scale shocks along the cavity walls. The SiO 
knot position is distant from the central protostars IRS 1 and IRS 
2 (see e.g. figure[^, which most likely prevents any UV heating 
from the central star to be operating. Furthermore, it lies out¬ 
side the envelope, as revealed by IRAC (|Tobin et al.|20101l, NH 3 
( 1 , 1 ) dBourke et all|1995a| fT997] l or N 2 H^ ( |Chen eFaLll^OOSl l 
emission. In other words. The SiO knot position offers the op¬ 
portunity to study pure small-scale shocks along cavity walls. 
However, the comparison of our shock models with those pre¬ 
sented in van Kempen et al. ( |2010b l or |Visser et al. ( |2012[ ) is not 
really meaningful. These authors indeed used a multiple-shock 
layers model, with a pre-shock density of 10 ^ ^ cm^^ close to the 
protostar, and 10 ^^ ^ cm“^ further away along the envelope (see 
figure 4 of|Visser et al.|201^. Furthermore the ID shocked lay- 
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Fig. 8. The CO excitation diagram produced from our best-fit model 
(red squares). Two temperature components can be fitted to our shock 
model over the PACS range of observations: a warm component (Tex 
= 216 K, dark blue line) fitting the level populations for Jap= 12 to 25 
(light blue squares), and a hot component (Tex = 422 K, light green line) 
fitting the level populations for Jup= 26 to 40 (dark green squares). 


ers they used were outputs of the Paris-Durham model we are 
also using in the present study. However they did not include 
the tip of the envelope that would correspond to the SiO knot 
position, as they solely focus on the outflow cavity in the vicin¬ 
ity of the protostar. The conclusion from our analysis is that our 
constrained pre-shock density of 10"^ cm~^ is consistent with the 
range of values they used (as they considered a decreasing den¬ 
sity further away from the protostar). 

Another, more classical approach to discuss energetics and 
physical conditions based on CO excitation consists of using 
excitation diagrams, as reviewed in |Visser| ( |2014| l. Figure 
presents the excitation diagram built from our best-fit model for 
CO. As extensively described in [Goldsmith & Langer| ( |1999| l, 
in local thermodynamical equilibrium conditions and under op¬ 
tically thin regime, the points in this type of diagram are ex¬ 
pected to fall on a straight line, whose inverse of the slope is the 
excitation temperature of the transitions (also see the descrip¬ 
tion of H 2 excitation diagrams in section 4.1 1 . As pointed out in 
jVisserj ( |2014| l, numerous studies can be found in the literature 
describing the building of this kind of excitation diagrams based 
on PACS observations centred on outflow-driving protostars of 
low- to high-masses, and their subsequent decomposition in at 
least two gas components, one warm (7’ex=320+50 K), and one 
hot (7’ex=820+150 K). Usually the breakpoint between these two 
components is around 1800 K, between the (25-24) and (26-25) 
transitions. We have produced a similar diagram from our best-fit 
model, as can be seen in figure]^ For comparison purposes, we 
have also applied a two-component linear fit of the lines accessi¬ 
ble to PACS, with the warm component (Tex = 216 K) fitting the 
level populations for 7up= 12 to 25, and a hot component (Tex = 
422 K) fitting the level populations for 7up= 26 to 40. The values 
we constrained for the excitation temperatures within these two 
components are systematically below those inferred from PACS 
observations centred on outflow-driving protostars of all possible 


mass (van Kempen et al. 2010a||Herczeg et al. 2012 

Goicoechea 

et al.|2012| Manoj et al.|2013[ Karska et al.|2013 

20T4r Green 

et al.||2013a|b[ Dionatos et al.||2013[ Lee et al.||2013[ Lindberg 


jet al.|2014| . This shows that pure shocks contribute to both the 
so-called ‘warm’ and ’hot’ components. Although the excitation 
temperature associated with the warm component is close to the 
observed values, the figure also shows that pure shocks fail to 
entirely account for any of these warm or hot components. In¬ 
deed the excitation temperatures that we find here are less than 
those derived from the observations in those papers. This means 
that close to the protostar, additional mechanisms, or different 
kinds of shocks should account for the CO observations. It also 
shows to which extent the collection of continuous temperature 
components that is generated by our shock models can be in¬ 
terpreted as a two-excitation component, through the analysis of 
CO excitation diagrams, for transitions with 7up=12 to 40. 

Outflow parameters. Finally, we studied the energetics asso¬ 
ciated with this pure shock position. Our first method is to follow 
the procedure presented in Anderl et al. ( 2014| l for the W44 SNR 
shock study. Indeed, we can infer the mass, momentum, and en¬ 
ergy associated with our best-fit model, under the assumption of 
a filling factor equal to one. The results are presented in the sec¬ 
ond column of table|^ First, the total mass contained in the beam 
is 1.8 X 10“^ Mq. This value is far greater than that determined 
for the mapped area (half of the inner part of the outflow) by N09 
(2.5xlO“^Mo), or Gial 1 (0.6x10“^Mq). This might be partly ex¬ 
plained by the fact that our beam size is rather large and that we 
decided to focus on the brightest H 2 spot in the map, contrary 
to those studies that operated on values averaged over larger ar¬ 
eas. More convincingly, this is probably due to the fact that our 
observations-models approach gives us access to the population 
of the (v = 0, 7 = 0, 1) H 2 levels, contrary to the observations. 
The contribution of these two levels to the total column density 
is significant: a linear fit to the Spitzer observations presented 
in Figure 1^ yields a column density A(H 2 )~ 6.9 x 10'® cm“^, 
whereas linearly fitting the (v = 0, 7 = 0, 1) part of the modelled 
rotational diagram yields A(H 2 )~ 1.8 x 10^* cm~^, which is 26 
times larger than that measured based only on the observations. 
Moreover, our value is consistent with the total mass of IMq 
for the northern lobe determined by jBourke et al.| ( [l997 1 based 
on CO (1-0) and (2-1) line observations (combined with *^CO 
and C'^O (1-0)). From a different perspective, the value that we 
found is simultaneously ~50 times smaller than that found in the 
bright CO positions of the W44 SNR studied in the same way 
by jAnderl et al.| ( [2014| l, where molecular shocks also propagate. 
The corresponding momentum is 0.4 Mq km s ', also a good 
factor ~100 below the values inferred using similar methods in 
W44. This is also a factor 10 below the value found based on 
more observational methods (described in the next paragraph) 
over a 12'.'5 beam in the massive star-forming region W28 A2 
that encompasses three different outflows at least (Gusdorf et al., 
in press). Finally, the dissipated energy is 4.2 x 10‘'^ erg, typi¬ 
cally two orders of magnitudes below the equivalent quantity in 
W44. This is also two orders of magnitudes below the energy 
dissipated in a 1275 beam in W28 A2 but comparable to what 
was found by [Gomez-Ruiz et all] ( |2013| l for the entire outflows 
(associated with low-mass YSOs) L1448-IRS3 and HH211-mm 
(based on the use of more observational methods). From this 
method, it seems that BHR71 can be defined as an energetic low- 
mass outflow, but the relatively high numbers we found could be 
the effect of the method we used, based on a sophisticated shock 
model. From the energetic point of view, the impact of the whole 
BHR71 outflow on the Galactic ISM is much less than that, e.g. 
of the whole W44 SNR, because the dissipated energy is smaller, 
but also because of its much smaller size: the BHR71 outflow is 
roughly covered by a rectangle of ~ 0.1 x 0.5 pc^, whereas a 
SNR such as W44 resembles a circle of ~26 pc radius. 
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We also used a second method, following [Bontemps et al.| 
( 1996| l; Beuther et al.| ( [200^ , to calculate the outflow parame¬ 
ters related to the considered species (CO, H 2 and SiO). These 
parameters (dynamical age fj, mass M, mass entrainment rate 
M, momentum P, mechanical force kinetic energy E^, and 
mechanical luminosity Lmech, are calculated using the equations: 




( 1 ) 


M - NX TtR^ X tMspecies, (2) 

M = MIta, (3) 


Table 3. Outflow parameters over the beam of our observations in H 2 , 
CO, and SiO. The 6vmax and values for H 2 were assumed to be identi¬ 
cal as for CO. The second column summarises the values extracted from 
the best-fit model (first method, see text). 


Species 

Model 

H 2 

CO 

SiO 

N (cm“^) 

3.7e21 

1.9e21 

3.2el7 

2.0el5 

M(Mo) 

1.8e-2 

1.3e-2 

1.5e-5 

1.5e-7 

dUmax (km s^') 

- 

40 

40 

40 

fd (yr) 

- 

1785 

1785 

1785 

M(Moyr-') 

- 

7.0e-6 

8.6e-9 

8.4e-ll 

P (Mq km s"') 

0.4 

5.0e-l 

6.1e-4 

6.0e-6 

Fm (Mq km s^' yr^') 

- 

2.8e-4 

3.4e-7 

3.4e-9 

Ek (erg) 

4.2e44 

2.0e44 

2.4e41 

2.4e39 

^mech (^0) 

- 

9.3e-l 

l.le-3 

l.le-5 


P - Mx6v„ 


(4) 


(5) 


= M X 5vl,J2, 


( 6 ) 


Tmech — 


(7) 


where N is the column density of the considered species, R the 
radius of our analysis region, and the approximate zero-base 
linewidth of the considered species. We used CO, H 2 , and SiO 
column densities extracted from our best-fit model associated 
with a filling factor of 1. The results are indicated in table 
Lacking the necessary spectral resolution for the H 2 data, we as¬ 
sumed that the 5urnax and fd values for H 2 are identical identical 
as for CO. This leads to overestimating the mass, momentum, 
and energy calculated based on the H 2 data with respect to the 
results purely obtained from the model (first method). This is be¬ 
cause the global shock velocity of our best-fit model, 22 km s“', 
is less than the zero base linewidth that we attributed to H 2 in 
this second method. Regardless, the CO mechanical luminosity 
we constrain is again consistent with that measured by |Bourke| 
|et al.| ( |1997| l for the whole outflow (0.5 Lq), whereas the corre¬ 
sponding value we calculated for H 2 exceeds that given by N09 
for the fraction of the outflow that they mapped (0.05 Lq)- 

We can now compare the outflow parameters as inferred 
from our observations of BHR71 with similar studies found in 
the literature for outflows associated with protostars of various 
masses. An impressive number of such studies have been aimed 
at isolated targets, making it difficult to give some perspective 


on our results based on a one-to-one comparison (e.g. Leurini 


e^D|20% Lebron et al. 2006jJFonteni_e^_aO2009J |Liu et alT 


2010( Guzman et al.|201~ [Cyganowski et al.|2011| l. We conse 

quently focus on a comparison with ‘survey studies’, i.e. studies 
that are aimed at extracting outflow parameters from a sample 
of sources. From this respect, we found that the outflow parame¬ 
ters derived from CO observations of BHR71 are obviously sys¬ 
tematically less than the values inferred in more massive envi¬ 
ronments (e.g. |Duarte-Cabral et al.|2013[ [Klaassen et al. 201 l|l. 


let alone in the outflows associated with O-type stars 


,opez- 


Sepulcre et al.|20()9ll. This is also true for SiO-related energetic 


parameters: BHR71 is less energetic from this point of view than 
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its massive counterparts ( |Sanchez-Monge et al.||2013] l. Indeed 
these authors computed the mass, momentum, and energy asso¬ 
ciated with the SiO (2-1) and (5^) lines. Their study was fo¬ 
cussed on the whole outflows, resulting in values considerably 
larger than those inferred here (4-110 and 2-35 Mq; 26-2130 
and II- 44 OM 0 km s^'; (0.2-75) and (0.06-16)xl04®ergs). On 
the contrary, the energetic parameters inferred from our study 
of one bright beam are similar to those inferred over the extent 
of the whole outflows associated with similar low-mass stars as 
IRS 1, for instance. This is the case for the outflows in the Perseus 
cloud (NGC1333, L1448 , HH211, L1455, e.g. [Curtis etal.|2010[ 
[Gomez-Ruiz et al.|2013| ). The conclusion is that either BHR71 is 
indeed more energetic than its low-mass counterparts, or that our 
method partly based on shock modelling naturally yields higher 
numbers than in all these studies, that often make use of fewer 
CO lines, for instance, than in the present work. 


6. Conclusions 

We have presented new observations of the BHR71 bipolar out¬ 
flows, obtained with SOFIA in '^CO (11-10), and (1^15), and 
with APEX in '^CO (3-2), (4-3), (6-5), (7-6), and ^^CO (3-2). 

We combined these data with existing datasets: in H 2 
(Spifzer/IRS) and SiO (APEX), and we have compared the ob¬ 
servations in the form of a CO flux diagram, an H 2 excitation 
diagram, and an SiO integrated intensity diagram to a grid of 
sophisticated shock models. 

Our best fit is non-stationary (C7-type) and has the following 
parameters: wh = 10^ cm“^, b - 1.5, - 22 km s“', and an age 

of 3800 years. The age and velocity of the shock model hint at 
the presence of more than one shock structure within the rather 
large beam of our observations, 24". This was also suggested 
by the fact that we had to assume different filling factor values 
for the different considered species. Erom the analysis of our 
ten best-fit models, we consider that the constrained values are 
quite robust, as these models all had hh = 10 ^ cm“^, b — I - 2, 
Us = 20 - 25 km s“', and an age of 3500^500 years. 

However this modelling can still be used to discuss the feed¬ 
backs of the shocks encompassed in our observations. We stud¬ 
ied its chemical feedback in terms of SiO chemistry, placing an 
upper limit of 4% of the total silicon-bearing material in the form 
of SiO in the grain mantles in the pre-shock region. 

We quantified the contribution of shocks to the excitation 
of CO around low-mass protostars surrounded by outflows, and 
shown that the CO excitation diagram from a shocked layer 
where the gas temperature is calculated point by point can be 
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interpreted as a two-component one for levels from 7i,p=12 to 
40. 

We inferred global outflow parameters from our shock 
model; a mass of 1.8 x 10 “^ Mq was measured in our beam, 
in which a momentum of 0.4 Mq km s“' was dissipated, corre¬ 
sponding to an energy of 4.2 x 10^^ erg. We also analysed the 
energetics of the outflow species by species. Both methods sug¬ 
gest that BHR7 lisa rather energetic outflow. 

Three perspectives lie ahead of the present study. The first 
is to generalise our analysis to the whole outflow, and to obser- 
vationally constrain the outflow energetic parameters over the 
whole mapped area. This will yield meaningful comparisons 
with observational studies of various similar objects. The sec¬ 
ond is to observe the SiO knot position with ALMA to resolve 
the multiple shock structures caught in our present single-dish 
beam. This should allow us to isolate a single shock structure, 
which would help us to get closer to the geometrical configura¬ 
tion that the Paris-Durham model was tailored to fit, and lead the 
shock analysis to a new level of understanding, both in terms of 
chemistry and energetics. Finally, at this point, more emission 
lines will also have to be included in our analysis, especially 
H 2 O lines. 
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Fig. A.l. The same field as in figure]^ shown in CO (6-5) (colours) 
convolved at the 24" resolution. The half-maximum contours of the 
emission of the CO (3-2) (black line), (4-3) (red), (6-5) (green), and 
(7-6) (blue) lines is also shown at the same 24" resolution. The SiO 
peak and knot, HH320 AB and knot 5 positions, as well as the beam 
sizes of the (6-5), (16-15) and (11-10) transitions are also indicated as 
in figure]^ 


Appendix A: Filling factors 

The half-maximum emission contours of the CO (6-5) and (3- 
2) lines were already shown in figure in thick black and red 
contours, at their nominal resolutions of 9" and 1871. For both 
the SiO peak and knot, the filling factor of the emission with 
respect to a beam of 24" (that of the CO (11-10) line observa¬ 
tions), inferred from the red, CO (3-2) contours is of the order 
of 1. As we performed our analysis over a 24" beam size, and 
not at the resolution of the CO (3-2) line, we decided to verify 
that this filling factor assumption was correct for all CO transi¬ 
tions at a 24" resolution. We hence convolved all the APEX data 
to this resolution. The result is shown in figure [AT| The dataset 
is remarkably consistent in terms of the size of the emitting re¬ 
gion: the half-maximum emission contour for all lines is broadly 
the same, except for the (7-6) transition, which could be due 
to insufficient signal-to-noise values. Based on this figure, we 
chose to operate with a filling factor of 1 for all CO transitions 
observed in the SiO knot position over a beam of 24". 


Appendix B: CO, SiO, and H 2 observables 

Appendix C: The CO flux diagram in the SiO peak 
position 


Table B.l. CO and SiO integrated intensity values, and H 2 level popu¬ 
lations measured over a beam of 24" centred on the SiO knot position. 
This means that the CO (16-15) (nominal resolution of 1673) and SiO 
(nominal resolution of 28") data were a posteriori corrected assuming 
the emission is extended over the largest beam. As such, the initial data 
were multiplied hy (16.3/24)^ and (28/24)^. The uncertainties are given. 
The integrated intensities were calculated between -4.5 and 50 km s“'. 
The CO and SiO integrated intensities correspond to a filling factor of 
1 . 
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Fig. C.l. The observational CO flux diagram over a beam of 24" ob¬ 
tained in the SiO peak position. The CO (16-15) line is not detected: an 
upper limit estimate is displayed in the form of a hlack arrow. 
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